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Abstract 
We show how accurate benchmark values of the surface formation energy of crystalline lithium hydride 
can be computed by the complementary techniques of quantum Monte Carlo (QMC) and wavefunction- 
based molecular quantum chemistry. To demonstrate the high accuracy of the QMC techniques, we present 
a detailed study of the energetics of the bulk LiH crystal, using both pseudopotential and all-electron 
approaches. We show that the equilibrium lattice parameter agrees with experiment to within 0.03 %, 
j3 , which is around the experimental uncertainty and the cohesive energy agrees to within around 10 meV per 

formula unit. QMC in periodic slab geometry is used to compute the formation energy of the LiH (001) 
surface, and we show that the value can be accurately converged with respect to slab thickness and other 
' O ' technical parameters. The quantum chemistry calculations build on the recently developed hierarchical 

C , scheme for computing the correlation energy of a crystal to high precision. We show that the hierarchical 

scheme allows the accurate calculation of the surface formation energy, and we present results that are well 
converged with respect to basis set and with respect to the level of correlation treatment. The QMC and 
hierarchical results for the surface formation energy agree to within about 1 %. 
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I. INTRODUCTION 

The surface formation energies of materials are key quantities in fields as diverse as nanotech- 
nology, mineral science and fracture mechanics. However, the accurate measurement of surface 
energies is fraught with difficulties, so there is often a need to rely on calculated values. In princi- 
ple, electronic-structure methods based on density-functional theory (DFT) should be capable of 
giving reliable surface energies, but in practice it is found that computed values depend strongly 
on the approximation used for the exchange-correlation energy.—""— There are two main kinds of 
electronic-structure technique that allow one to go beyond DFT and achieve better accuracy: quan- 
tum Monte Carlo (QMC), and the wavefunction-based correlation techniques usually associated 
with molecular quantum chemistry (QC). We show here how these two very different approaches 
can be used in a complementary way to produce accurate benchmark values for surface formation 
energy, using as a test case the (001) surface of crystalline LiH. 

There have been many DFT calculations of the surface formation energies a of different kinds 
of materials, including ionic compounds, covalent semiconductors and metals. In some cases, the 
variation of predicted a values with the assumed exchange-correlation functional has been studied, 
and it is found that generalized gradient approximations (GGA) such as PBE and PW91 often give 
a values that are ~ 30 % lower than those predicted by the local density approximation (LDA).— ^— 
Since GGAs are generally more accurate than LDA for bonding energies, and since the energy 
needed to form a surface would seem to be closely related to the energy needed to break bonds, 
it might be expected that GGA values of a would be more accurate. However, in the few cases 
where there are reliable experimental data, this expectation is not fulfilled, and the rather scattered 
evidence suggests that the LDA may be more accurate.— "^i^ A connection has been made with the 
superiority of LDA over GGAs for the surface energy of jellium.— 

In this rather confused situation, it is helpful to seek ways of computing benchmark values 
of a which do not suffer from the uncertainties of DFT. Quantum Monte Carlo, and specifically 
diffusion Monte Carlo (DMC) '^ offers one way of achieving this. It is well established that 
DMC is usually much more accurate than DFT for the energetics of extended systems, and 
there are ways of systematically improving its accuracy. Nevertheless, it is subject to errors that 
are not completely controllable, and this is where the methods of molecular quantum chemistry 
can play an important role. The electron-correlation techniques that we use here are mainly 
second-order M0ller-Plesset theory (MP2) and the coupled-cluster scheme CCSD(T) (including 
single, double and a perturbative treatment of triple excitations). Efforts to apply these QC 
techniques to the energetics of extended systems go back many years, particularly using the 
so-called incremental approacb^"^. More recently, the MP2 approximation has been implemented 
for periodic systems in several code o^^'^^ . The present authors have reported a technique 
referred to as the hierarchical method for applying molecular QC methods to perfect crystals, 
and for the case of LiH have shown that it can deliver the cohesive energy to an absolute ac- 
curacy of ~ 5 meV per formula unit and the equilibrium lattice parameter to better than 0.1 %.— »^ 

We have chosen to study the surface energetics of LiH, partly because it is a material for which 
we expect DMC to give very high accuracy, and partly because we already know that hierarchical 
QC is very accurate.— The crystal has the rock-salt structure, and the simplicity of this structure 
facilitates the calculations. We have several main aims. First, we want to show that DMC does 
indeed deliver high accuracy for the properties of the LiH crystal, particularly if we use all-electron 
rather than pseudopotential DMC (pp-DMC). Second, we report our periodic slab calculations of 
a for the LiH (001) surface, using both pseudopotential and all-electron DMC, and we show that 
we can achieve a high degree of convergence with respect to slab thickness and other technical 
parameters. Third, we show that the hierarchical QC scheme that gives such good accuracy for 



bulk LiH also provides a practical way of obtaining benchmark values of a. The hierarchical 
methods allow us to calculate explicitly the contribution of core-valence correlation to a, and 
we shall see that this is significant. Naturally, close agreement between the a values computed 
by the QMC and QC approaches supports the credibility of both, and this will be carefully assessed. 



II. TECHNIQUES 

A. Quantum Monte Carlo 

For present purposes, the name quantum Monte Carlo refers to two techniques for determining 
the ground-state energy of a many-electron system (for reviews, see e.g. Refs.— i^). Our high- 
precision results are obtained using diffusion Monte Carlo (DMC), a technique that projects out 
the ground state by evolving the many-electron wavefunction in imaginary time with the aid of 
an approximate trial wavefunction. An optimized form of this trial function is computed using 
variational Monte Carlo (VMC), which is an implementation of the variational principle of quantum 
mechanics. The VMC and DMC calculations in this work are performed using the CASINO 
package^. 

The trial wavefunctions used here have the standard single-determinant Slater-Jastrow form: 

*T = D^D^e-^ , (1) 

where D^ and D\^ are up- and down-spin Slater determinants of single-electron orbitals V'n- Electron 
correlation is approximately described by J, which is a sum of three types of terms: electron- 
electron terms u, electron- nucleus terms Xi ^iid electron-electron- nucleus terms /. These three 
terms contain parameters that are optimized using VMC, so as to make ^t as close as possible to 
the true ground-state wavefunction. The optimization works with the local energy Ei^ = ^^ H"^^, 
where H is the many-electron Hamiltonian. We follow the common procedure of varying the 
parameters so as to minimize the variance of E]^ (the variance would be zero if \I't were the exact 
ground-state wavefunction). VMC can be used equally well as an all-electron technique or with 
non-local pseudopotentials to represent the interaction between valence electrons and atomic cores. 

The idea of DMC—'^'iS, is to represent the exact many-electron wavefunction <i> as a density 
of brownian particles, or 'walkers'. In the evolution of the wavefunction according to the time- 
dependent Schrodinger equation in imaginary time, the optimized approximation \I't from VMC is 
used to guide the walkers, in a manner related to importance sampling. DMC aims to stochastically 
simulate the diffusion, birth, death and drift of the walkers, which, after an equilibration period, 
samples the exact ground-state wavefunction. In practice, the fermionic nature of electrons prevents 
DMC from being completely exact, and the nodal surfaces of the wavefunction are constrained 
to be those of ^t ^ this is the well-known fixed-node approximation^^. We shall see that this 
approximation incurs only small errors in the present work. A number of other technical issues have 
to be addressed, including time-step errors, pseudopotential errors, the choice and representation 
of the single-electron orbitals V'n, and the stability of walker populations, and we summarize these 
next. The treatment of system size errors will be discussed in Sec. IIIBi 

The walkers propagate by using the approximate small-time-step Green's function as a transi- 
tion probability in configuration space. The approximate Green's function also includes a term that 
gives a probability for a given walker to 'branch' (become two walkers) or to be discarded entirely. 
The use of a discrete time-step incurs errors, but these can be rendered negligible by the usual 
procedure of extrapolating to the zero-time-step limit. We shall present both pseudopotential and 
all-electron DMC calculations on the LiH bulk and surface. For the pseudopotential work, we use 



the Dirac-Fock pseudopotentials due to Trail and Need o^^'^^ . It is difficult to treat non-local pseu- 
dopotentials in DMC, and we employ the usual locality approximation^^, which introduces errors 
proportional to the square of the difference between ^t and the exact ground-state wavefunction 
$. The comparison of our pseudopotential and all-electron results will help us to quantify these 
errors. 

The single-electron orbitals ipn used in the trial wavefunction $t (see Eq. ([T])) were generated 
by DFT calculations with the LDA functional. We make this choice because there is considerable 
evidence^^i^ that this gives a ^t that is closer to the true ground state. The ^pn were computed 
by plane-wave calculations with the quantum espresso package^^. However, the direct use of ipn 
in a plane- wave representation in DMC is very inefficient, and instead we re-expand the ipn in a 
blip- function (B-spline) basis^'^, using the standard relation between the blip-grid spacing and the 
plane- wave cut-off. In the case of all-electron DMC, a further modification is necessary, since it 
is crucially important that ^t has the correct electron-nuclear cusp at the nuclear positions. The 
technique we have used to ensure this with the blip basis is described in Appendix A. 

Since walkers can branch or be discarded after each step, the walker population fluctuates. A 
reference energy in the approximate Green's function allows us to bias the branching, and thus 
control the population. However, in regions of particularly low energy (especially divergences 
at point charges), this mechanism is not enough, and a walker trapped in this region (and its 
offspring) can branch repeatedly, causing a population explosion which destroys the statistics of 
subsequent moves. 



B. QMC for bulk and surface energies 

Correction for errors due to the limited size of the periodically repeated cell is important in the 
calculation of both bulk and surface energies. As usual, we distinguish between single-particle and 
many-body errors. The former are due to the fact that /c-point sampling cannot be performed with 
DMC, and are analogous to those that would arise in single-particle methods such as DFT without 
/c-point samplng; the latter are due to the spurious interaction of electrons with their periodic 
images. To correct for single-particle errors, we use the formula-: 

where E^^ and E^^^ are the energies of the given cell with DMC and DFT (no /c-point sampling 
with DFT), E^'^ is the DFT energy of the cell with perfect /c-point sampling, and E^o is the 
corrected DMC energy; if enough data for different system sizes are available, a can be treated as 
a fitting parameter. 

One way of correcting for many-body size errors is to use a modified form of the Coulomb 
interaction known as the Model Periodic Coulomb (MFC) interaction in the DMC calculations.—"— 
We used this technique, in combination with Eq. ([2]) for our all-electron calculations on bulk LiH. 
An alternative approach is the scheme due to Kwee et al.—, which corrects for both single-particle 
and many-body errors in a single formula: 

rp _ pDMC I pLDA pKZK /o\ 

which somewhat resembles Eq. ([2]). Here, E^^ is a DFT-like energy of the cell (no A;-point 
sampling), which uses a functional designed to mimic the sum of single-particle and many-body 
errors, while E^^ is the same as E^"^ in Eq. ([2]), evaluated with the LDA functional. We used 
this scheme of Kwee et al. for the pseudopotential calculations on the bulk. Whichever method 



is used to correct for the many-body size errors, in the case of the bulk calculations we apply a 
further two-point extrapolation to remove residual finite-size errors. This extrapolation employs 
the formula: 

E^ = {NEm-MEm)/{N-M) , (4) 

where Ei\f and Em are the DMC energies per formula unit of supercells containing N and M 
formula units respectively.— i^^i"— 

Our DMC calculations of the surface formation energy are performed in slab geometry, so that 
we work with slabs having infinite extent in the plane of the surface and having a specified number 
N of ionic layers. Periodic boundary conditions are applied in the surface plane, so that we have 
supercell geometry only in two dimensions. With the blip basis set used for the present work, it is 
unnecessary to apply periodic boundary conditions in the direction perpendicular to the surface. 

As usual, the surface formation energy a is the work needed to create an area A of new surface, 
starting from the perfect bulk crystal, divided by A. In slab geometry, if £'siab(-^) is the energy per 
supercell of the A^-layer slab, Vs\ah{N) is the number of formula units per supercell of the A^-layer 
slab, and Cbuik is the energy per formula unit of the bulk crystal, then a is given by: 

a = lim (^siab(A^) - '^siab(A^)ebuik) M > (5) 

where A is the total surface area (both faces) per supercell of the slab. In Eq. ([5]), we must take the 
limit as the number of layers A^ and the surface dimensions of the supercell both tend to infinity. 
For comparison with experimental data, the ionic positions in the slab should also be relaxed to 
equilibrium, but in the present work we are concerned mainly with comparing different theoretical 
approaches, and we focus on the unrelaxed value of cr, for which all ions in the slab have their bulk 
positions. 

Instead of using Eq. ([5j) directly, we prefer to use the well-known procedure of extracting a 
from a series of slab calculations of increasing N , using the fact that as A^ — t- c«, -Esiab has the 
asymptotic form: 

E,uh{N)^Aa + NEuyer- (6) 

Here, a is the surface formation energy with the chosen surface supercell, and £^iayer is the bulk 
energy per ionic layer with this supercell. Eq. Q is equivalent to Eq. ([5]), but the extraction of a 
for a given surface supercell from Eq. ([6]) is usually more robust. We note that the value of -Eiayer 
can be cross-checked against independent calculations on the bulk crystal, since A^-E'iayer/i^siab(A^) 
should be very close to Cbuik- 

When correcting for finite-size errors in slab geometry, compensation for many-body errors 
poses technical problems, and we therefore used only the single-particle correction of Eq. ([2]). 

C. Correlated quantum chemistry 

We show here how the hierarchical metho d^^'^^ , originally developed to treat bulk crystals, can 
be used to calculate surface formation energy. We recall that the hierarchical method begins by 
separating the total energy e*°* per primitive cell of a crystal into Hartree-Fock and correlation 
parts: 

gtot ^ gHF ^ gcorr _ ^7) 

The correlation energy e™^^ is further separated into a molecular contribution and the so-called 
"correlation residual" : 



^corr „corr 



e3 + Ae™-. (8) 



In the case of a compound AB having the rock-salt structure, e™^^ is the correlation contribution to 
the binding energy of the AB molecule, with the bond length taken equal to the nearest-neighbour 
distance in the crystal. 

The hierarchical method works by combining energies of a sequence of finite cluster o^^'^^ in 
such a way as to eliminate surface effects. For the rock-salt structure, we take cuboidal clusters 
having /, m and n ions along the three perpendicular edges. By conventional quantum chemistry 
techniques, we can compute accurately the total energy E\^^ of each I x m x n cluster, which is 
then decomposed into Hartree-Fock, molecular and residual parts: 

i??°*„ = eZ. + limn eZ\ + ^EfZ ■ (9) 

The total energy per primitive cell in the infinite crystal is then: 

e*°*= lim J_£;[°t„ = eH^ + en + , Hm j^AEfZ ■ (10) 

l,m,n^oo Lmn l,m,n^oo tmn 

We calculate the Hartree-Fock contribution e^^ using standard periodic codes, and e'^l\ is obtained 
by conventional quantum chemistry techniques. To perform the limiting process in the third term 
on the right, the hierarchical method expresses AEf^ as: 

AEfZ = 8E°°° + 4[(/-2) + (m-2) + (n-2)]SO°i 

+ 2[(m - 2)(n - 2) + (n - 2)(/ - 2) + (/ - 2)(m - 2)]^;°^^ 

+ {l-2){m-2){n-2)E^^^ , (11) 

where the coefficients E^^^, E^^^, E^^^ and E^^^ represent the energies of corner, edge, face and 
bulk sites, respectively. [Note that the definitions of E^^^, E^^^ and E^^^ are affected by our 
decision to use factors {I — 2){m — 2)(n — 2), (m — 2)(n — 2), etc, rather than lmn, nin, etc. The 
reason for making this particular choice of factors is discussed in Ref.— .] 

Our procedure for obtaining the values of the coefficients in the limit of infinite l,ni,n requires us 
to extract E^^^, E^^^ , E^^^ and E^^^ from sets of four independent clusters, and then systematically 
to increase the size of the clusters in these sets, as described in detail in Ref. Il3l . For the cohesive 
energy, only the limiting value of E^^^ is needed, since Ae'^™'^ = 2E^^^. However, the procedure 
also yields the limiting values of E^^^, E^^^ and E^^^. The value of E^^^ can be used to obtain the 
value of the unrelaxed surface formation energy a. 

The coefficient E^^^ is the contribution to the energy of a large cluster from an atom in the 
interior; E^^^ is the same for an atom on the surface. When a surface is formed by opening a 
gap in the crystal, each atom in the newly formed surface contributes E^^^ — E^^^ to the energy 
difference. The area of the surface occupied by each atom is a^/2, so the correlation contribution 
to the formation energy of a new surface is 2{E^^^ — E^^^)/a'^ per ion. We therefore obtain 

a = a«^ + 4(E0ii-i^"i)/a2. (12) 

We compute the Hartree-Fock part a^^ using standard periodic codes (details will be given later). 

D. Zero-point corrections 

Both quantum Monte Carlo and quantum chemistry techniques employ static calculations, and 
ignore zero-point energies. In this work, in order to facilitate comparison with experiment, all bulk 
calculations are corrected for zero-point energy. These corrections are calculated using DFT and 
the linear response method. The PBE functional^ is used, since this is known to give accurate 
phonon frequencies; our tests with the LDA functional showed little change in the zero-point 
energy. The calculations were performed using the Quantum Espresso package^^. 



III. BULK LIH WITH QUANTUM MONTE CARLO 

We present first our pseudopotential DMC calculations on the bulk, which already give quite 
high accuracy, and also provide valuable information about the effect of system-size errors on the 
cohesive energy E^oh (the energy per formula unit relative to free atoms), the equilibrium lattice 
parameter gq and the bulk modulus B. The all-electron DMC bulk calculations reported at the 
end of this Section will show that an explicit treatment of core-valence correlation improves the 
accuracy still further. 



A. Pseudopotential calculations 

The Dirac-Fock non-local pseudopotentials^^'^* that we use are rather hard, and we found that 
a plane- wave cut-off of 4080 eV and a correspondingly fine blip-grid spacing was needed to produce 
accurate orbitals. It proved straightforward to eliminate DMC time-step errors: a time-step of 
0.025 au reduced the error below 1.0 meV per formula unit, which is much greater accuracy than 
we need. 

To study system-size errors, we calculated the DMC total energy for several values of the atomic 
volume, using cubic supercells containing 3x3x3 and 4x4x4 primitive crystal cells (54 and 
128 atoms). In addition, a DMC calculation on the 5x5x5 system (250 atoms) was performed 
at a single atomic volume. The correction of Kwee et al— was then applied to each calculation 
using Eq. ^, and two-point extrapolation (Eq. ^) was used to reduce the remaining many-body 
finite-size errors. To illustrate the effect of system-size errors, we show in Fig. [T] plots of ii^coh as a 
function of atomic volume from our DMC calculations on the 3x3x3 and 4x4x4 supercells, as 
well as the results corrected for size errors, compared with our earlier quantum-chemistry cohesive 
energies obtained with and without core-valence correlation effectsi^. In each case, we show also 
the third-order Birch-Murnaghan fit to the results, 

E{V) = i^o - ^ ((4 - B',)^^ - (14 - 3i?^)^ + (16 - 3i^^)^) , (13) 

where £"0, Vq and Bq represent the equilibrium energy, volume and bulk modulus, respectively, and 
Bq is the first derivative of the bulk modulus. 

The values of Ecoh, oo and Bq obtained from the fit are given in Table HI Several important 
points emerge from these results. First, the system-size effects consist almost entirely of a vertical 
shift, i.e. a constant energy offset, of the -E'coh(^) curves, so that they cause only small errors in oq 
and Bq. For example going from the 3x3x3 supercell to the extrapolated system changes qq by 
only ~ 0.1 %. Second, comparing the raw values of Ecoh from DMC on the 3x3x3 and 4x4x4 
supercells to the fully-corrected and extrapolated value suffers from substantial errors of 133 and 
52 meV, the KZK correction reduced these to 50 and 22 meV respectively. Third, the DMC value 
of -Ecoh) even after correction, still disagrees with the quantum chemistry value of -Ecoh without core 
correlation energy by ~ 36 meV. This last point indicates that the effect of the pseudopotential 
approximation must be significant. In order to make further progress, all-electron DMC is needed, 
and we report on this next. 

B. All-electron calculations 

In order to perform accurate all-electron DMC on the LiH crystal, we have to address several 
technical challenges. First, as noted in Sec. Ill A[ the trial wavefunction must accurately satisfy the 



FIG. 1: Cohesive energy as a function of primitive cell volume for the extrapolated pp-DMC calcualtions and 
the quantum chemistry calculations with and without core-valence correlation effects. DMC calculations on 
finite supercells are included to show convergence. The lines indicate a Birch-Murnaghan fit to the data. 
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TABLE I: Calculated bulk properties with both pseudopotential and all-electron DMC and hierarchical 
quantum chemistry with and without core effects. The cohesive energy is calculated at 4.084A in all cases. 
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/ k Bo / GPa ^coh / eV 
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DMC Extrap. 
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"Including the KZK correction 

'This value is extrapolated to infinite-size, zero-timestep using six separate calculations with 3x3x3 and 4x4x4 



supercells and timesteps of 0.004, 0.002 and 0.001 a.u. 



Kato cusp condition at the nucleus, in order to ensure stability of the walker population. Second, 
because of the rapid variation of the orbitals near the nucleus, extremely fine blip-grids are needed. 
Third, we expect to need much shorter time-steps than for pseudopotential calculations. 

The technique used to ensure that the cusp condition is satisfied is outlined in Appendix A. 
One symptom of the rapid variation of the orbitals near the nucleus is the slow convergence of 
the DFT total energy with respect to plane-wave cut-off in the quantum ESPRESSO calculations 
used to generate the orbitals. Given the high memory requirements caused by the high cut-offs 
in the DMC calculations, we took the orbitals to be sufficiently converged when they produced 
stable walker populations. For our final all-electron DMC calculations, we used orbitals generated 
using the LDA functional, with a plane-wave cut-off of 6.8 x 10^ eV. The associated blip-grid had 
a spacing of half the natural grid dictated by the plane-wave cut-off. 

We made detailed tests on the time-step needed to ensure the accuracy of the all-electron 
calculations. In Fig. [21 we show the results of tests on the 3x3x3 supercell, showing how the 
total energy converges with respect to time step. The figure also shows a linear fit to the results, 
which is clearly adequate, as expected from earlier work.~ A time step of 0.004 au gives an error 
of only 10 meV/formula unit, which more than suffices to give accurate results for the equilibrium 
ao and Bq. 



FIG. 2: Total energy vs. time-step. This shows the convergence of the ae-DMC calculations with respect 
to time-step and is for a 3x3x3 supercell. A linear extrapolation to zero timestep is included. 
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In order to obtain the best possible results for the cohesive energy as a function of volume 
Ecoh{V), we used the fact made clear in Sec. IIIIAI that results with the 3x3x3 supercell differ 
only by an almost constant energy offset from results converged with respect to supercell size, 
this offset being in the region of 30 meV. Our procedure in the all-electron DMC calculations was 
therefore to calculate EcohiV) first with the 3x3x3 supercell and a time step of 0.004 au. We 
then added a constant correction energy to these results, obtained from DMC calculations with 
time steps of 0.001, 0.002 and 0.004 au, all performed on both the 3x3x3 and 4x4x4 supercells. 
At each time step, the usual two-point extrapolation (Eq. (jl])) to infinite supercell size was made, 
and a final linear time-step extrapolation was then made. 

The cohesive energy curve from all-electron DMC is compared in Fig. [3] with our pseudopo- 
tential DMC curve and the results from quantum chemistry. The resulting equilibrium values of 
-E^coh) Qq and Bq are compared in Table HI We see that the all-electron DMC value of oq agrees 
with the experimental and quantum chemistry values to within ~ 10^^ A (0.03 %). Values of Bq 
are much more difficult to obtain accurately, but the all-electron DMC value agrees with quantum 
chemistry to within ~ 4 %, and both are reasonably consistent with experimental values, which 
span a range of ~ 15 %. The all-electron DMC and quantum chemistry values of the equilibrium 
ii^coh differ by 13 meV/formula unit. The quantum chemistry value is believed to be somewhat 
more accurate than this, so that some of this 13 meV may be due to fixed-node error. 



IV. SURFACE FORMATION ENERGY OF LIH WITH QMC 



The methods of Sees. IIIAI and IIIBI have been used to calculate the formation energy of the 
LiH (001) surface, first with pseudopotentials, then with all-electron DMC. All the calculations 
were done with the lattice parameter oq = 4.084 A. 



A. Pseudopotential calculations 



Exactly the same pseudopotential methods were used for the calculations on slabs as were used 
in the bulk calculations of Sec. IIIIAi and the trial orbitals were generated using DFT with the 
LDA functional, as before. These orbitals were then re-expanded in B-splines using a spacing 
corresponding to T^/2kma.x where A;max is the modulus of the largest plane- wave vector. For each 



FIG. 3: Cohesive energy vs. primitive cell volume. Here the DMC energy differences have been calcu- 
lated using a 3x3x3 supercell and timestep of 0.004 a.u. and the point at 17.03A'^ calculated using the 
extrapolation procedure in the text. 
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surface supercell and each number N of ionic layers, the Jastrow factor of 1-, 2- and 3-body terms 
was re-optimized using variance minimization. 

To extract the values of a and i^iaycr for each chosen surface unit cell, we performed calculations 
of the total slab energy -E'siab(-^) foi' numbers of ionic layers from 3 to 6 using a 4 x 4 surface 
unit cell (18 ions per layer in the repeating supercell). Single-particle size errors were corrected for 
using Eq. ([2]) with a set equal to 1. Table HIl shows the convergence of a with respect to the slabs 
used when fitting to Eq. ([6]). We have also performed calculations for slabs 3 and 4 using a 3 x 3 
surface unit cells (18 ions per layer in the repeating supercell). Comparing directly the £73^4 from 
the two different surface unit cells differed by only 0.006(5) Jm~^ indicating the finite-size error. 
The resulting best value for a from the pseudopotentlal calculations is 0.373 Jm~^. 

TABLE II: pp-DMC surface formation energy calculated using slabs of different thicknesses. Calculations 
performed on 4 x 4 surface unit cells. 



Slabs used ct / Jm 



3,4,5,6 

4,5,6 

5,6 



0.369(2) 
0.373(3) 
0.379(6) 



B. All-electron calculations 

The all-electron DMC techniques used for the slab calculations were essentially the same as those 
used for the bulk (Sec. lIIIB]) . However, the memory requirements for the B-spline coefficients were 
so much greater than for the bulk that we had to reduce the plane-wave cut-off used to generate 
the orbitals from 6.8 x 10^ to 3.4 x 10^ eV. This primarily made the DMC runs more susceptible 
to population control issues and resulted in a higher statistical error on the final values compared 
to the pseudopotentlal work. For the same reason, were able to perform all-electron calculations 
only for the 3x3 surface unit cell, and the largest number of ionic layers that we could handle was 
N = 5. We know from the pseudopotentlal calculations that the finite size errors are under control 
using 3x3 surface unit cells assuming the LDA correction is used. The introduction of the tightly 
bound core electrons is not expected to increase the finite-size errors. 
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A further technical issue in the all-electron slab calculations was concerned with the optimization 
of the Jastrow factor. In order to obtain wavefunctions that produced stable DMC runs we found it 
necessary to optimize the Jastrow factor using the energy minimisation scheme within VMC. This 
tended to increase the variance of the local energy of the trial wavefunction slightly with respect 
to variance minimisation. However it did reduce the number of population explosions during the 
DMC runs. 

The timestep adopted (0.004 a.u.) was the same as used in the all-electron bulk work, since 
the tests done there indicated that this is sufficient. Our final all-electron DMC result for the 
surface formation energy is a = 0.44(1) Jm.~^. We note that the explicit inclusion of Li core states 
increases o" by ~ 0.07 J m~^, which is a significant effect at the level of accuracy sought in this work. 



V. HIERARCHICAL QUANTUM CHEMISTRY FOR SURFACE ENERGY 

The formation energy of the of LiH (001) surface was also computed using quantum chemistry 
techniques. The Hartree-Fock component u^^ was determined from slab calculations, see Eq. ([6]), 
using the CRYSTAL and VASP codes. The effect of electron correlation was accounted for using 
the hierarchical method as described in Section [II CI 

Both CRYSTAL and VASP employ periodic boundary conditions, so that the calculations are 
performed on an infinite array of slabs, with a vacuum gap separating successive slabs. The vacuum 
gap was chosen to be 26 A, large enough to ensure that there is no interaction between neighbour- 
ing slabs. Careful attention was also paid to convergence with respect to fc-point sampling and 
basis-set completeness. A previous high accuracy Hartree-Fock study of bulk LiH was performed 
using CRYSTAL by Paier et al^. The basis set described in that work was used for the present 
calculations. In order to ensure basis set completeness, layers of "ghost" atoms were added above 
and below each surface. The "ghost" atoms were basis functions centred on the sites of atoms in 
the next layer but without the nuclei or electrons. The convergence of a^^ with respect to these 
"ghost" atoms was tested using slabs of four and five layers, see Table IIIII The introduction of 
the "ghost" atoms has a significant effect on a^^ and two layers are necessary to achieve basis set 
completeness; this number of layers was used in all our calculations. 

TABLE III: Convergence of cthf using CRYSTAL with respect to the number of layers of "ghost" atoms 
above and below the surface. Based on two-point extrapolations from slabs of 4-5 layers. All energies are 
quoted in Jm~^. 



Ghost layers cthf 






0.43835 


1 


0.19886 


2 


0.19849 


3 


0.19854 



Calculations on slabs of two to eight layers were performed using both periodic codes, and the 
method outlined in Sec. IIV Al was used to extract values of o"^^. The resulting values are shown 
in Table |Vl We note that the VASP value is slightly lower than the CRYSTAL value. Since 
CRYSTAL provides a direct all-electron calculation, and we have established that the CRYSTAL 
result is converged with respect to basis set, we suggest that the VASP value may be a slight 
underestimate. This may be due to the PAW potentials used: the standard PBE potentials were 
used, and while harder potentials are available for H, it was not possible to reach convergence with 
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these potentials. Previous studies of bulk LiH with VASP have reported small discrepancies in 
the Hartree-Fock result^. The CRYSTAL results converge with respect to slab thickness to give 
a value of 0.198(1) Jhi'^. 

TABLE IV: Hartree-Fock approximation to surface formation energy for LiH, a = 4.084 A, using CRYSTAL 
and VASP. All energies are quoted in Jni^^. 



Slabs used CRYSTAL VASP 



2-8 


0.20005 


0.19363 


3-8 


0.19883 


0.19114 


4-8 


0.19825 


0.19001 


5-8 


0.19819 


0.18944 


6-8 


0.19836 


0.18926 


7-8 


0.19703 


0.18864 



The correlation component of the surface formation energy was calculated using the hierarchical 
method. The convergence of the hierarchical coefficients is shown in Fig. |4] and using the methods 
described in Ref. [ij the values can be converged to within a few tenths of niEh. In brief, a 
reference calculation was performed using A^ = 64 and frozen-core MP2 theory in the cc-pVTZ 
basis set. Corrections for core correlation (Jcore), basis-set incompleteness (5basis) and higher- level 
correlation treatments ((5CCSD(T), 5CCSDT, 5CCSDT(Q)) and the diagonal Born-Oppenheimer 
correction were also computed using smaller basis sets and smaller values of N. The use of smaller 
values of A^ for small corrections was validated in earlier work^^ . Complete details of the hierarchical 
results are given in Table |Vl An error of 0.2 mEh in the hierarchical coefficients corresponds to 
an error of 0.005 Jm~^ in the surface formation energy. To facilitate comparison with both DMC 
results, a^"" has been calculated with and without correlating the core electrons. 



FIG. 4: The convergence oi E'^^^ (red circles) and E^^^ (blue squares) with respect to maximum cluster size 
N using MP2/cc-pVTZ for LiH a = 4.084 A 

-3.5 
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VI. DISCUSSION 



We summarize in Table IVII our QMC and hierarchical quantum chemistry results for the for- 
mation energy a of the LiH (001) surface, and we compare with the predictions of DFT using the 
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TABLE V: Correlation contributions to the surface formation energy of LiH at a = 4.084 A and total 
calculated surface formation energies with and without core correlation. In each row, the N value specifies 
the maximum number of ions in the hierarchical calculation. Correlation-consistent basis sets have been 
used throughout, and cc-p(C)VXZ is abbreviated (C)VXZ. The CCSDT and CCSDT(Q) calculations were 
performed using MRC C'^^'^^ as a module in Molpro and the diagonal Born-Oppenheimer correction (DBOC) 
calculations were performed using PSI3^^. All other calculations were performed using Molprc'^*. 

reference E°^^/mEh g^Vm^h cr'^°"7Jm-^ details 

+0.1886 MP2/VTZ iV = 64 

+0.0321 MP2/CVTZ - MP2/VTZ iV = 16 

+0.0000 MP2/V[T,Q]Z iV = 36 

+0.0004 MP2/V[Q,5]Z iV = 16 

+0.0091 CCSD(T)/VTZ - MP2/VTZ iV = 16 

+0.0063 CCSDT/VDZ - CCSD(T)/VDZ A^ = 8 

+0.0008 CCSDT(Q)/VDZ - CCSDT/VDZ A^ = 8 

-0.0015 HF/VTZ A^ = 8 





-4.196 


-6.000 


Score 


-0.840 


-1.147 


5basis 


+0.268 


+0.268 




+0.162 


+0.158 


(5CCSD(T) 


+0.625 


+0.538 


(5CCSDT 


-0.149 


-0.209 


(5CCSDT(Q) 


-0.015 


-0.023 


DBOC 


+0.041 


+0.055 



frozen core 
"total 



0.2037 Sum terms above except Score 
0.2358 Sum all terms above 



,HF 



—Static 

frozen core 

static 
"total 



0.198 

0.402 
0.434 



LDA, PBE and rPBE functionals, these DFT results being taken from our earlier work—. The very 
close agreement between the DMC and QC results for a confirms that both approaches give high 
accuracy, and shows that the results can be used as benchmarks for assessing DFT approximations. 
The DFT values of a span a remarkably wide range, with LDA overestimating it by 7 %, and PBE 
and rPBE underestimating it by 23 % and nearly 40 % respectively. It is interesting to note that 
the Hartree-Fock value a^^ of 0.20 J m^^ (Table |V]) accounts for less than half the full value of 
a, so that the importance of an accurate treatment of correlation is clear. It is also noteworthy 
that valence-core correlation gives a surprisingly significant contribution of -^ 10 % to a. We noted 
in the Introduction the scattered evidence that LDA tends to give better values of a than GGA 
approximations, and the present work shows that this is the case for LiH. 

TABLE VI: Calculated surface formation energy with both pseudopotential and all-electron DMC and 
hierarchical quantum chemistry with and without core effects. The calculations are performed at 4.084A 
in both cases. DFT data from previous work is included. The DFT values are for the lattice parameters 
optimized with the given functionals.'* 



method 



a I Jm 



DMC pseudopotential 0.373(3) 

DMC all-electron 0.44(1) 

Quantum Chemistry (froz. core) 0.402 

Quantum Chemistry (with core) 0.434 

DFT LDA 0.466 

DFT PBE 0.337 

DFT rPBE 0.272 



An important part of the evidence that our DMC and hierarchical QC calculations give results 
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of benchmark quality for a is the very high accuracy of the two completely independent approaches 
for the energetics of the LiH crystal; for hierarchical quantum chemistry, this was already shown 
in detail in Ref. 13, and a substantial part of the present paper has been devoted to showing 
the same thing for QMC. In both approaches, the calculated cohesive energy is correct to ~ 
15 meV/formula unit, and the equilibrium lattice parameter to within better than 0.1 %. It is 
clear from both sets of calculations that an adequate treatment of core-valence correlation must 
be included. The reliability of the calculated values for a is confirmed by the very close agreement 
(within ~ 0.01 Jm^^) between the values given by the two approaches. Here too, core- valence 
correlation is important, giving a contribution of ~ 0.03 Jm^^ to a. 

In both theoretical approaches, a key technical issue is system size effects. In QMC, the calcu- 
lations are done directly in periodic boundary conditions, but large repeated systems are needed 
because we rely on F-point sampling. For the bulk crystal, we have shown that significant correc- 
tions need to be made for size errors, but that these are successful in reducing the errors in the 
cohesive energy to ~ 15 meV/formula unit. For the DMC slab calculations of a, we have presented 
evidence that the calculated a is very well converged (to within 0.006 Jm~^) with respect to both 
slab thickness and size of the surface repeating unit. In the hierarchical approach, the Hartree-Fock 
part of a is calculated in periodic boundary conditions, and we have shown that size errors in the 
slab calculations can be made negligible. Size errors in the correlation residual contributions are 
well controlled in the hierarchical scheme, and, as can be seen from Fig. [H are of the order of 0.1 
m£Jh per ion. 

It would now be timely to extend the present calculations to other materials. In fact, we 
have reported QMC calculations of a for MgO (001) several years agc^, though it might be worth 
repeating the calculations with the improved pseudopotentials now available. Our hierarchical 
quantum chemistry scheme can be applied without change to MgO and other materials having the 
rock-salt structure, and we hope to report both QMC and hierarchical calculations on LiF in the 
near future. However, it is important to emphasise that the hierarchical scheme also works well for 
materials having other crystal structures, so that there is now rather wide scope for using it, with 
or without QMC, for calculations of a. We remark that for some materials it will be essential to 
include the effects of surface relaxation. This is not a significant issue for most rock-salt materials, 
and we have shown^ that for LiH relaxation reduces a by only ~ 0.01 Jm~^. But in other cases 
(corundum is a famous example), relaxation makes a large difference to a, and would probably 
need to be estimated from DFT calculations. 

To conclude, we have shown that: (a) quantum Monte Carlo calculations give extremely 
accurate results for the energetics of the LiH crystal, particular when Li core electrons are 
explicitly included, and there is excellent agreement with results from the hierarchical quantum 
chemistry scheme; (b) these two independent techniques give almost identical benchmark results 
for the formation energy of the LiH (001) surface; (c) the benchmark value of a lies between DFT 
predictions from the LDA and GGA approximations, the LDA value being somewhat better than 
GGA. 
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Appendix A: Generalised cusp correction 



A scheme for modifying real orbitals expanded in a Gaussian basis set so that they satisfy the 
Kato cusp conditions^^''^^ is described in Ref. |4l|. In the present work we make use of an extension 
of this scheme which allows the Kato cusp conditions to be imposed on complex orbitals expanded 
in any smooth basis set. 

For simplicity, we restrict our attention to the case of an all-electron nucleus of charge Z at the 
origin in the following discussion. In the scheme of Ref. l4ll . the s-type Gaussian basis functions 
are replaced by radial functions in the vicinity of the nucleus. These functions impose the cusp 
conditions and make the single-particle local energy resemble an "ideal" curve that was found 
empirically to be satisfied by a wide range of Hartree-Fock atomic orbitals. In the scheme used 
in this work, instead of replacing part of the orbital, we add a spherically symmetric function of 
constant phase to the orbital. The function added to uncorrected orbital V'(r) is 



A-(/'(r) = exp(i0o) (t){r) — (f){r) 6(rc — r), 
where Q is the Heaviside function, Tc is a cutoff length, Oq = arg['i/;(0)], 

"exp(-i6'o) 



(r) = Re 



47r 



sphere 



V'(r) dQ. 



and 



= C + exp(ao + air + a2r + a^r + a^r 



(Al) 



(A2) 



(A3) 



where C is a real constant and the {a} are real constants to be determined. In practice (j){r) 
is evaluated by cubic spline interpolation; the spherical averaging of the uncorrected orbital is 
performed on a radial grid at the outset of the calculation. C is chosen so that (j){r) — C is positive 
everywhere within the Bohr radius of the nucleus. 
The uncorrected orbital may be written as 



■0(r) = exp(i6lo)0(r) + r?(r). 



(A4) 



where r/(r) consists of the / > spherical harmonic components of ip{r), together with the phase- 
dependence of the I = component. Note that exp(i^o)i;^(0) = V'(O)) and hence 7y(0) = 0. We may 
now apply the scheme of Ref. |4l| to determine {a} and r^ with exp(i0o)'?^ a-nd exp(i0o)'^ playing 
the roles of the uncorrected and corrected s-type Gaussian functions centered on the nucleus at 
the origin. The constant phase exp(i6'o) cancels out of the equations that determine the {a} [Eqs. 
(9)-(13) in Ref. |41], so the determination of the {a} and Tc is exactly as described in Ref. l4ll . 
except that we do not need to modify Z when more than one nucleus is present, because r/(0) = 0. 
Suppose the orbital is of Bloch form ip{r) = Mk(r)exp(ik • r), where Uk has the periodicity 
of the primitive cell. Let {Rp} be the set of primitive-cell lattice points. The orbital may be 
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corrected at each all-electron nucleus in the primitive cell at Rp = using the scheme described 
above. The phase of the orbital (and hence cusp-correction function) at the corresponding nucleus 
in the primitive cell at Rp 7^ is exp(ik • Rp) times that for the primitive cell at Rp = 0; the 
cusp-correction function is otherwise identical. The corrected orbital is of Bloch form. 
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